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1.  INTRODUCTION 
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Collision  theory  is  well  established  as  an  elementary  classical  theory 
for  bimolecular  gas  reactions,  and  the  concepts  are  readily  incorporated  in 
the  direct  simulation  method.  It  is  a phenomenological  approach  in  that  the 
binary  reaction  rate  is  obtained  as  the  product  of  the  kinetic  theory 
collision  rate  for  collisions  with  energy  in  excess  of  the  activation 
energy  and  the  probability  of  reaction  or  steric  factor  P . The  steric 
factor  may  be  regarded  as  the  ratio  of  the  reactive  cross-section  to  the 
total  collision  cross-section  and  is  inferred  from  a comparison  of  the 
theoretical  prediction  with  the  measured  reaction  rate.  The  simplest 
theory  (for  example*)  employs  a constant  steric  factor  and  considers  only 
the  translational  collision  energy  in  relation  to  the  activation  energy. 
Refinements  of  the  method  employ  a collision  energy  dependent  steric  factor 
and  allow  for  the  possible  contribution  of  the  internal  energy  of  the 
molecules  in  the  collision.  Collision  theory  may  also  be  extended  to  ter- 
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molecular  reactions  (for  example  ) by  assigning  a lifetime  to  each  binary 
collision  and  regarding  the  triple  or  ternary  collision  as  a further  binary 
collision  between  the  pair  of  molecules  in  the  initial  collision  and  a 
third  molecule. 

The  reaction  of  most  obvious  interest  in  an  aerodynamic  context  is 
the  dissociation-recombination  of  oxygen  and  nitrogen.  These  are  bi- 
molecular reactions  with  finite  activation  energies  for  dissociation  and 

termol ecu! ar  reactions  for  recombination.  The  direct-simulation  method  was 
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applied  to  this  reaction  by  Bird.  These  calculations  employed  a constant 
steric  factor  and  the  results  for  the  dissociation  rate  were  in  agreement 
with  the  theoretical  prediction.  Also,  the  density  and  temperature  depen- 
dence of  the  equilibrium  state  was  in  agreement  with  the  law  of  mass 
action  for  this  reaction.^  The  gas  model  was,  however,  unrepresentative 
of  a real  diatomic  gas  in  that  there  was  no  internal  energy.  A further 
difficulty  was  that  computing  time  considerations  ruled  out  calculations 
at  temperatures  small  compared  with  the  characteristic  dissociation  tem- 
perature. Similar  calculations  which  include  the  effects  of  internal  energy 
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modes  have  recently  been  made  by  Larsen. 
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Bird*’  developed  a collision  theory  analysis  of  the  dissociation- 
recombination  reaction  to  produce  an  expression  for  the  equilibrium 
state  that  is  identical  in  form  to  the  law  of  mass  action.  If  a parti- 
cular diatomic  gas  is  to  be  simulated,  a comprison  of  the  theoretical 
and  experimental  dissociation  rates  gives  the  appropriate  steric  factor, 
while  a comparison  of  the  theoretical  and  experimental  values  for  the 
equilibrium  degree  of  dissociation  provides  the  binary  collision  lifetime 
parameter.  The  model  corresponding  to  this  extended  theory  was  combined 
with  Borngakke  and  Larsen's  model ^ for  the  internal  energy  modes  and  was 

used  to  compute  the  dissociation  relaxation  zone  behind  a strong  shock 
8 

wave  in  nitrogen.  The  ability  to  calculate  the  equilibrium  state  was 
particularly  valuable  in  setting  up  the  boundary  conditions  for  the  shock 
wave  simulation.  The  simulation  results  were  in  good  agreement  with 
corresponding  experimental  results. 

A difficulty  that  was  encountered  in  the  nitrogen  simulation  was 

that  the  negative  temperature  exponent  in  the  rate  equation  was  below  the 

lower  limit  allowed  in  the  theory.  This  had  to  be  allowed  for  by  a 

temperature,  as  well  as  collision  energy,  dependence  of  the  steric  factor. 

The  following  section  presents  a further  extension  of  the  bimolecular 

theory  to  include  the  contribution  of  the  internal  energy  modes.  As  well 

as  being  more  general,  it  allows  much  larger  negative  temperature  exponents 

in  the  rate  equations  and  would  have  coped  with  the  accepted  value  for 

nitrogen.  Unlike  the  existing  collision  theory  based  on  the  "effective 

1 3 

square  term"  approach  , it  is  not  restricted  to  temperatures  that  are 
small  in  comparison  with  the  characteristic  temperature  of  the  reaction. 

A further  extension  Is  that  the  model  now  deals  with  inverse  power  law 
molecules  rather  than  being  restricted  to  hard  sphere  molecules.  The 
wider  applicability  of  the  new  model  is  particularly  useful  when  dealing 
with  cases  which  involve  a very  large  number  of  reactions. 
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2.  TECHNICAL  APPROACH 

2.1  DERIVATION  OF  STERIC  FACTORS  FROM  REACTION  RATES 
A typical  bimolecular  reaction  may  be  written 

A + BXC  + D , (1) 

where  A,  B,  C,  and  0 represent  separate  molecular  species.  The  rate 
equation  for  species  A may  be  written 

- = kf(T)  nAng  - kr(T)  ncnQ  . (2) 

where  n is  the  number  density  with  the  subscripts  representing  the  species, 
t is  the  time  and  k^(T)  and  kr(T)  are  the  forward  and  reverse  rate  co- 
efficients. The  temperature  T dependence  of  the  rate  coefficient  is 
generally  of  the  form 

k(T)  = a Tb  exp  (-Ea/kT)  , (3) 

where  a and  b are  constants,  Ea  is  the  activation  energy,  and  k is  the 
Boltzmann  constant.  The  problem  is  to  find  an  expression  for  the  steric 
factor,  or  collision  probability,  Pr  as  a function  of  the  collision  energy 
Ec  that  leads  explicitly  to  Eq.  (3)  in  the  context  of  collision  theory. 

The  collision  energy  Ec  is  equal  to  the  sum  of  the  relative  transla- 
tional energy  in  the  collision  and  the  fraction  internal  energy  of  the  two 
molecules  that  may  contribute  to  the  reaction.  The  translational  energy 
is  equal  to  1/2  mrcr  where  mr  = m^/firij  + mg)  is  the  reduced  mass  and 
cr  is  the  relative  speed.  The  internal  energy  is  best  dealt  with  by 
specifying  that  the  number  of  degrees  of  freedom  of  the  two  molecules  that 
may  contribute  energy  to  the  reaction  is  ^ and  Cg.  This  number  may  range 
from  zero  to  the  number  of  degrees  of  the  molecule,  although  there  are 
restrictions  based  on  angular  momentum  considerations  on  the  number  of 

3 


FT 


rotational  degrees  of  freedom  that  can  be  active  in  a reaction.  It  is 
readily  shown  that,  for  an  equilibrium  gas,  the  distribution  function  of 
Ec/kT  is 
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where  n is  the  exponent  of  inverse  power  law  intermolecular  force  and  r 
is  the  Gamma  function.  Collision  theory  assumes  that  the  forward  rate 
equation  is  given  by 
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where  is  the  equilibrium  collision  rate  for  a species  A molecule  with 
species  B molecules,  and  E^  is  the  activation  energy  such  that  the  reaction 
probability  is  zero  for  E_  less  than  E . Now 
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where  o-j.  is  the  total  collision  cross-section  and  the  bar  denotes  the  mean 
value.  A comparison  of  Eqs.  (2)  and  (5)  then  shows  that 
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For  hard  sphere  (n  = »)  molecules,  the  standard  result  is  that 
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where  d is  the  molecular  diameter  while,  for  general  values  of  n, 
^ - Wj  (K/m^  , 


(8) 


where  WQm  is  the  cut-off  value  of  the  dimensionless  impact  parameter  and 
K is  the  constant  in  the  force  law.  The  parameter  Wgm  is  related  to  the 
lower  cut-off  value  that  must  be  chosen  for  the  deflection  angle  in  the 
collision.  The  value  of  Pr  will  depend  on  this  choice  and  the  generalized 
collision  theory  is  unsatisfactory  when  the  overall  aim  is  merely  the 
calculation  of  the  steric  factor  P . However,  this  is  not  a difficulty 
when  this  factor  is  used  in  the  direct  simulation  method  since  a cut-off 
value  must  be  specified  anyway. 

If  a constant  steric  factor  is  employed  in  Eq.  (6),  it  may  be  taken 
outside  the  integration  which  then  yields 
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For  integer  values  of  the  first  argument,  the  incomplete  Gamma  function 

may  be  written  as  a power  series  in  E /kT.  This  is  effectively  done  in 

a 1 

the  "activation  in  many  degrees  of  freedom"  theory  and  the  "square  term" 
theory.  All  but  the  leading  term  of  the  power  series  may  be  discarded 
when  kT  « E.  but,  while  this  is  satisfactory  for  the  dissociation- 

a 

recombination  reaction,  it  is  out  of  the  question  for  fast  exothermic 
reactions.  This  problem  may  be  avoided  by  the  introduction  of  a steric 
factor  proportional  to  some  power  of  E -E„  multiplied  by  the  term 

C a 


This  term  enables  the  incomplete  Gamma  function  to  be  transformed  to  a 
Gamma  function  that  does  not  involve  the  temperature,  and  leads  to  the 
desired  form  of  Eq.  (3)  for  k(T).  Moreover,  this  modification  has  physical 
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as  well  as  mathematical  advantages.  It  means  that  Pr  need  no  longer  be 
constant  for  very  large  values  of  E and  the  rise  from  zero  at  E_  = E, 
is  no  longer  discontinuous.  Similar  forms  for  the  steric  factor  have, 
in  fact,  been  proposed  for  these  reasons. 

The  full  expression  for  the  steric  factor  that  yields  Eq.  (3)  when 
substituted  into  Eq.  (6)  with  Eqs.  (4)  and  (8)  is 


(9) 
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The  symmetry  factor  e is  equal  to  1 for  unlike  molecules  and  2 when  both 
collision  partners  are  of  the  same  species.  The  use  of  Eq.  (7)  in  place 
of  Eq.  (8)  gives  the  following  simplified  form  of  Eq.  (9)  for  hard  sphere 
molecules: 


77T7T2-  Vnk  / 
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The  minimum  value  of  b is  now 


!) 


rather  than  the  - allowed 


by  earlier  theory  which  does  not  take  into  account  the  contribution  of 
the  internal  energy  modes.  The  extended  theory  would  have  covered  the 
experimental  values  for  nitrogen  that  caused  difficulties  in  the  shock 

O 

wave  calculations.  This  general  behavior  is  consistent  with  the  more 
restricted  "square  term"  theory,  in  that  increasing  negative  value  of  b 
are  associated  with  a greater  contribution  from  the  internal  degrees  of 
freedom. 
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The  use  of  Eqs.  (9)  or  (10)  to  obtain  the  steric  factor  for  use  in 
a simulation  program  ensures  that  the  reaction  rates  will  agree  with  those 
from  the  continuum  rate  equation  in  flow  situations  where  near  transla- 
tional and  rotational  equilibrium  prevail  and  both  the  discrete  particle 
and  the  continuum  approaches  are  applicable.  The  discrete  particle 
approach,  through  the  direct-simulation  method,  permits  results  to  be 
obtained  for  chemically  reacting  flow  situations  with  marked  non-equilibrium 
in  the  translational  and  internal  modes.  The  critical  question  is  the 
degree  of  realism  in  the  modeling  of  the  real  flow  in  this  regime.  It 
should  be  kept  in  mind  that  the  above  equations  for  the  collision  proba- 
bility constitute  phenomenological  models  which  are  not  unique  and  which 
have  been  based  largely  on  considerations  of  mathematical  tractability. 

The  ideal  would  be  to  have  tabulations  of  reaction  cross-sections  as  func- 
tions of  the  energy  states  and  impact  parameters  based  on  quantum 
calculations  and  experiment.  The  progress  in  this  direction  has  been 
reviewed  by  Toennies,1*  and  very  little  of  this  information  is  currently 
available  for  reactions  of  engineering  interest.  When  a classical  approach 
based  primarily  on  continuum  rate  data  must  be  used,  the  above  model 
appears  to  offer  the  most  practical  approach  for  chemically  reacting  flows 
in  the  transition  regime.  Hopefully,  the  physical  basis  is  sufficiently 
sound  for  the  predictions  to  be  at  least  qualitatively  valid. 

In  order  to  consider  the  direct  simulation  method  as  an  alternative 
to  conventional  continuum  calculations,  the  following  conditions  should  be 
assessed. 

The  advantages  are: 

• The  gas  need  not  be  in  an  equilibrium  or  quasi-equil ibrium 
state.  Not  only  may  each  species  have  a separate  temperature, 
but  the  velocity  distributions  functions  may  vary  drastically 
from  the  Maxwellian  distribution  and  the  internal  modes  may 

be  out  of  equilibrium  with  the  translational  mode. 

• The  calculations  make  direct  use  of  the  reactive  cross- 
sections  as  functions  of  the  molecular  properties  and  collision 
impact  parameters.  The  fact  that  these  do  not  have  to  be 
averaged  to  give  rate  constants  as  a function  of  temperature 

is  particularly  advantageous  for  electron-neutral  and  electron- 
ion  reactions  for  which  the  cross-section  versus  energy  curves 
often  have  closely  spaced  peaks  and  valleys.  The  vibrational 
states  may  be  regarded  as  reactions.  Radiative  reactions  may 
also  be  included. 
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The  disadvantages  are: 
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• The  number  of  representative  simulated  modules  is  generally 
restricted  to  some  thousands  and  the  total  number  of  simulated 
collisions  to  some  millions.  This  means  that  it  is  impractical 
to  simulate  reactions  which  occur  less  frequently  than  once 
every  thousand  elastic  collisons  (i.e.,  steric  factors  less 
than  about  10~3).  This  rules  out  some  very  slow  bimolecular 
reactions  and  generally  rules  out  termolecular  reactions. 

• There  are  sampling  difficulties  associated  with  trace  species, 
although  these  may  be  alleviated,  to  some  extent,  through  the 
use  of  weighting  factors  such  that  the  number  of  real  molecules 
represented  by  each  simulated  molecule  varies  from  species  to 
species. 

• Most  of  the  currently  available  data  is  presented  as  rate 
coefficients  as  a function  of  temperature,  rather  than  as 
reactive  cross-sections  as  a function  of  the  collision  energy. 
This  data  must  be  analytically  converted  to  give  the  same 
answer  in  the  quasi-equilibrium  situations  for  which  the  con- 
tinuum approach  is  valid.  (The  magnitude  of  the  steric  factors 
given  by  this  conversion  process  provides  a useful  check  on 
the  physical  validity  of  the  rate  data.) 

The  above  considerations  lead  to  the  conclusion  that  direct  molecular 
simulation  is  most  suitable  for  rapidly  varying  transient  phenomena  in 


the  upper  regions  of  the  atmosphere.  It  can  be  used  in  situations  where 
the  conventional  continuum  approach  would  be  invalid. 

Illustrative  Example 

The  direct  simulation  method  has  been  applied  to  the  vibrational 
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excitation  of  steam.  The  reason  for  choosing  steam  is  that  it  has  an 
anomalously  rapid  vibrational  relaxation  rate  for  V-T  collisions.  The 
scatter  in  the  data  is  such  that  the  temperature  exponent  is  uncertain 
and  a constant  value  of  3.5  x 10"17m3s'1per  molecule  has  been  used  for 
the  deactivation  rate.  The  characteristic  temperature  of  the  reaction 
is  2293K  so  that  the  activation  rate  coefficient  given  by  the  product  of 
the  deactivation  coefficient  and  exp  (-2295/T).  The  resulting  steric 
factors  from  Eq.  (10)  are  shown  in  Figure  1. 

A homogeneous  relaxation  case  was  computed  for  initial  conditions 
of  T = 650K,  n = 1023  nf 3 and  n*  = 0.  Here,  n and  n*  are  the  number  den- 
sities of  the  ground  state  and  first  vibrational  level,  respectively. 

Eq.  (2)  becomes 
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and  energy  conservation  requires, 


T = 650  - 918  (n*/n)  . (12) 

Figure  2 compares  the  simulation  results  with  a numerical  solution  of 
eqs.  (11)  and  (12).  The  two  are  in  good  agreement  even  though  this  is  a 
sufficiently  fast  reaction  for  a substantial  departure  from  local  transla- 
tional and  rotational  equilibrium  to  be  expected.  The  equilibrium  concen- 
tration ratio  is  also  in  good  agreement  with  the  equilibrium  value  even 
though  the  simulation  procedures  based  on  collision  theory  cannot  be 
expected  to  satisfy  detailed  balancing. 

In  order  to  illustrate  the  effects  of  translational  equilibrium  in  a 
transition  regime  flow,  the  steam  model  was  then  applied  to  the  steady 
flow  through  a screen  or  grid.  This  flow  is  schematically  represented  in 
Figure  3 which  also  defines  the  flowfield  dimensions  and  the  initial  and 
boundary  conditions.  This  is  a one-dimensional  flow  in  which  the  grid  is 
a hypothetical  device  which  reflects  a specified  percentage  of  the  mole- 
cules passing  through  it.  The  initial  temperature  was  chosen  such  that 
the  initial  equilibrium  excitation  is  negligible  and  the  grid  solidity 
is  sufficiently  large  to  produce  a temperature  with  a significant  degree 
of  equilibrium  excitation.  The  temperature  profiles  along  this  flow  ahead 
of  the  grid  shown  in  Figure  4 are  very  similar  to  those  for  a normal  shock 
wave.  There  is  a temperature  jump  at  the  grid  location  and,  since  the 
downstream  Mach  number  in  this  flow  is  barely  supersonic,  there  is  some 
upstream  influence  from  the  open  end.  The  parallel  temperature  Tx  based 
on  the  velocity  components  in  the  flow  direction  overshoots  the  overall 
temperature  T by  a significant  margin.  The  normal  temperature  Tn  based  on 
the  other  components  overshoots  T by  a small  margin,  while  the  rotational 
temperature  Tr  lags  behind  T. 

Figure  5 compares  the  simulation  results  for  the  density  ratio 
n*/(n  + n*)  with  those  obtained  from  the  solution  of  eq.  (11)  using  the 
values  of  stream  velocity,  density  and  overall  temperature  given  by  the 
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simulation.  The  non-equilibrium  efffects  are  significant  and,  as  would  be 
expected  from  the  temperature  profiles  in  Fig.  4,  they  are  greater  when 
the  collision  energy  E is  obtained  from  the  translational  modes  only. 

Other  calculations  showed  that  the  vibrational  overshoot  is  much  larger 
if  the  rate  coefficient  is  proportional  to  a finite  positive  power  of 
temperature,  but  is  not  greatly  affected  by  the  replacement  of  the  basic 
hard  sphere  model  by  an  inverse  ninth  power  law  model. 

The  above  example  illustrates  the  scope  of  the  direct-simulation 
method  when  chemical  data  is  available  only  in  the  form  of  bulk  reaction 
rates.  Approximate  quantum  cross-section  predictions  are  available  for 
T-V  reactions  such  as  that  for  steam  in  the  second  example.  The  direct- 
simulation  approach  would  also  be  particularly  useful  for  the  reverse  task 
of  deducing  rate  coefficients  or  even  flowfields  from  these  predictions 
in  order  to  make  comparisons  with  experiment. 

2.2  NEGATIVE  ION  COMPOSITION  IN  THE  IONOSPHERE 

A set  of  chemical  reactions  representative  of  the  important  processes 
in  the  ionosphere  has  been  investigated  by  the  direct  simulation  Monte 
Carlo  method  described  earlier.  In  this  calculation,  a homogeneous  gas 
mixture  initially  containing  electrons,  atomic  oxygen,  molecular  oxygen, 
and  molecular  nitrogen.  Equation  (10)  is  used  to  convert  the  temperature 
dependent  rate  coefficients  given  below  to  collision  energy  dependent 
steric  factors. 

k(cm'Vs) 

< 10"16  for  T < 2000°K 

(1.4±2)x10-«(300)exp(^0)(ffli) 

for  (195  < T < 600) 

Exchange:  0"  + 02  -*•  0 + 0^  < 10"10 

0~2  + 0 - o2  + o'  ?io-10 


Formation:  e + 0g  -*■  0"  + 0 

e + 02  + 02  -*■  0~2  + 02 
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Detachment:  0"  + 0 -*•  02  + e 


O'  + 02  - 0 + 02  + e 


1.7  + 0.3  * 10~10 


2.3  x 10'9  exp  (^20)  for  T.  < 

i . / i 


02  + 0 -*•  Og  + e 


02  + 02  - 02  + 02  + e 


02  + N2  - 02  + N2  + e 


0‘  + N2  -*■  0 + N2  + e 


3.2  x 10’ 


(2.71-3)  MO’10  (4)  exp  (4»0) 


for  (375  < T < 600) 


(1.941.4)  xjo-12^)1'5^  (4p) 
for  (375  < T < 600) 

3.4  x io-9  exp  (=2600) 


Most  of  the  rate  coefficients  were  obtained  from  the  DNA  Reaction  Rate 
Handbook  (DNA. 19484),  March  1972. 

Two  formation  reactions  are  considered.  One  of  these 
e + 02  + 02  -*■  02  + 02 

is  a termolecular  reaction  with  a rate  coefficient  that  indicates  that 

it  occurs  only  once  every  107  binary  0,,  09  collisions  at  a number  density 
q in  ^ ^ J 

of  10  to  10  cm  . The  other  reaction 


e + 02  •+  0"  + 0 

is  bimolecular  but,  being  dissociative,  has  a high  activation  energy  of 
approximately  3.6  eV.  Above  this  energy,  the  steric  factor  is  of  the 
order  of  0.2.  The  exchange  and  detachment  reactions  have  steric  factors, 
in  the  energy  range  of  interest,  that  range  from  the  order  of  10"3  to 
unity. 


The  test  case  that  has  been  chosen  is  an  atmosphere  with  a number 
of  10“  cm”3  (approx.  200  km),  with  initial  mole  fractions  of  0.08  electrons, 
0.74  atomic  oxygen,  0.14  molecular  oxygen  and  0.04  molecular  nitrogen, 
and  at  a temperature  of  2000°K.  The  energy  of  the  electrons  is  instan- 


taneously raised  to  3.6  eV  and  the  reactions  are  followed  for  0.4  seconds. 

The  composition  is  shown  in  Figure  6 as  a function  of  time.  The 
mole  fraction  of  0"  increases  to  almost  0.04  around  0.12  seconds  and  then 
decreases.  The  fraction  of  0^  continues  to  increase  over  the  period  of 
the  calculation.  Figure  7 shows  the  species  temperature  over  the  same 
period.  They  are  still  a way  from  equilibrium  at  the  end  of  the  calculation 
(but  this  is  a comparatively  short  calculation).  In  Figure  8 the 
differences  in  the  translational  and  internal  temperatures  shown  as  a 
function  of  time  indicate  the  nonequilibrium  effects  of  the  internal 
degrees  of  freedom. 
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3.  CONCLUSION 


The  application  of  the  direct  simulation  Monte  Carlo  method  to 
describe  the  chemistry  in  the  ionosphere  is  potentially  a very  useful 
computational  tool.  Although  the  calculations  presented  are  restricted 
to  a reaction  cell  simulation  where  there  are  no  spatial  gradients  in 
the  thermodynamic  variables  the  method  does  describe  the  effects  of  the 
electrons  and  ions  and  nonequilibrium  influences  which  are  not  accounted 
for  in  the  continuum  descriptions.  The  procedure  can  be  used  to  para- 
metrically study  the  effect  of  high  energy  electrons  on  the  chemistry 
of  the  ambient  ionosphere  or  as  a screening  procedure  for  possible 
reaction  chains.  This  computational  procedure  will  be  incorporated 
in  the  future  into  the  general  Monte  Carlo  ion  collection  code  developed 
under  previous  contracts.  The  resulting  code  will  then  be  capable  of 
describing  the  disturbance  of  the  ionosphere  due  to  the  passage  of  a 
charged  rocket  borne  probe,  including  the  effects  of  molecular  collisions, 
electric  field  effects,  and  nonequilibrium  chemistry. 
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Figure  3.  A SCHEMATIC  REPRESENTATION  AND  DATA  DEFINITION 
OF  THE  TEST  CASE  FOR  THE  FLOW  OF  STEAM  THROUGH 
A GRID 
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Figure  4.  TEMPERATURE  PROFILES  FOR  THE  FLOW  OF  STEAM  THROUGH  A GRID 
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